cap log close
log using "$results/3c_additional_analysis_tables.log", replace
	
********************************************
**# Fig 1: Import heatmap & initial tariffs
********************************************	
qui {

**# Panel A: Import heatmap	

**## Setup
import excel "raw/sic87_3digit.xlsx", sheet("sic2") clear first
drop if missing(sic2)
save "$processing/sic2_w_desc.dta", replace
	
use "analysis/SIC87_regression_Dataset.dta", clear

keep if year==1979

global tariff_controls  other_tariff_change  dln_pstar_vw_7279_alt  mfa rho
global prod_controls  	ln_cap_lab_78 skill_emp_share_78 women_share 
global rti_controls  	ln_invest_def_78 lag_invest_def_ch rti   automation78   mat_ship_def_78

egen full = rowmiss($tariff_controls $prod_controls $rti_controls ave_ols_ch9 ave_iv_ch9  )
keep if full==0

keep sic
merge 1:m sic using "raw/xm_sic87_72_105_20120424.dta",
drop if wbcode ==""
drop if sic==.
gen sic2 = int(sic/100)
rename customs imports
collapse   (sum) imports, by(wbcode year sic2)

preserve 
	keep if year==1978
	collapse (sum) imports,  by(wbcode)
	gsort  -imp
	gen c_rank = _n
	replace c_rank = 11 if c_rank>10
	replace wb = "OTH" if c_rank>10		
	collapse c_rank (sum) imports,  by(wbcode)
	egen timp = total(imports)
	gen c_share = round(imports/timp,.005)
	format imports % 6.0fc
	drop timp
	tempfile top_countries
	save `top_countries', replace
restore
	
preserve
	keep if year==1978
	collapse  imports  , by(sic2)
	gsort  -imp
	gen s_rank = _n
	replace s_rank = 11 if s_rank>10
	replace sic2 = 100 if s_rank>10
	collapse s_rank  (sum)  imports, by(sic2)
	
	egen timp = total(imports)
	gen s_share = round(imports/timp,.005)
	format imports % 6.0fc
	drop timp
	tempfile top_sic
	save `top_sic', replace
restore

keep if inlist(year, 1979, 1988)
merge m:1 wb  using `top_countries' , keepus(wbcode) keep(1 3)
replace wb = "OTH" if _m==1
drop _m
collapse  (sum) imports, by(sic2 wb year)

merge m:1 wb  using `top_countries' 
drop _m
assert c_share!=.
sort c_rank

merge m:1 sic2 using `top_sic' , keepus(sic2)  keep(1 3)
replace sic =100 if _merge==1
drop _m

merge m:1 sic2 using `top_sic' 
drop _m
collapse c_rank s_rank c_share s_share (sum) imports, by(sic2 wb year)
	
merge m:1 sic2 using  "$processing/sic2_w_desc.dta", keep(1 3) nogen
replace desc = "Other" if desc == ""	

egen id = group(wbc sic)
xtset id year, delta(9)
gen ln_imports = ln(imports)
gen dln_imports = round(f.ln_imports - ln_imports,.05)
keep if year==1979
sort c_rank s_rank

gen rank = . 
loc cutoffs .25 .5 .75 1 1.15 1.25 1.5  2 4 10 
loc rank = 1
foreach c in `cutoffs' {
	replace rank = `rank' if dln_imports<`c' & rank==.
	loc ++ rank	
}

sort c_rank 
gen nc_rank = 12 - c_rank
drop c_rank
rename nc_rank c_rank

format dln % 03.2fc
gen x_pos = 0.0
gen y_pos = 0.0
gen x_pos2 = 0.1
gen cy_pos = c_rank-.5
gen s_lab_pos = s_rank + 0.25
	
gen s_sh = round(s_share*100,0)
tostring s_sh, replace force
replace s_sh  = " (" + s_sh + "%)"
replace desc = desc + s_sh
gen c_sh = c_share*100
tostring c_sh, replace force
replace c_sh = "(" + c_sh + "%)"
	
	
**## Make graph
	
replace rank = 0 if dln_<0

tw scatter c_rank s_rank if rank==0, m(S) msize(huge) mlc(black) mc(blue%40) ///
|| scatter c_rank s_rank if rank==1, m(S) msize(huge) mlc(black) mc(blue%20) ///
|| scatter c_rank s_rank if rank==2, m(S) msize(huge) mlc(black) mc(blue%10) /// 
|| scatter c_rank s_rank if rank==3, m(S) msize(huge) mlc(black) mc(blue%05) /// 
|| scatter c_rank s_rank if rank==4, m(S) msize(huge) mlc(black) mc(purple%05) /// 
|| scatter c_rank s_rank if rank==5, m(S) msize(huge) mlc(black) mc(purple%10) /// 
|| scatter c_rank s_rank if rank==6, m(S) msize(huge) mlc(black) mc(purple%20) /// 
|| scatter c_rank s_rank if rank==7, m(S) msize(huge) mlc(black) mc(red%40) /// 
|| scatter c_rank s_rank if rank==8, m(S) msize(huge) mlc(black) mc(red%50) /// 
|| scatter c_rank s_rank if rank==9, m(S) msize(huge) mlc(black) mc(red%60) /// 
|| scatter c_rank s_rank if rank==10, m(S) msize(huge) mlc(black) mc(red%70)  /// 
|| scatter c_rank x_pos  if s_rank==1  , m(i) mlab(wb) mlabangle(0) mlabpos(3) mlabc(black) mlabsize(vsmall) ///
|| scatter cy_pos x_pos2  if s_rank==1  , m(i) mlab(c_sh) mlabangle(0) mlabpos(3) mlabc(gs2%40) mlabsize(vsmall) ///
|| scatter y_pos s_lab_pos  if c_rank==11 , m(i) mlab(desc) mlabangle(90) mlabpos(3) mlabc(black) mlabsize(small) ///
ylabel(none) xlabel(none) ///
ytitle("") xtitle("") ///
caption("(a) 1979-1988 Import Growth by Country and Sector", position(6) justification(center) size(medium)) ///
yscale(lcolor(gs2%0)) ///
xscale(lcolor(gs2%0)) ///
legend(order(1 "<0.00" 2 "<0.25" 3 "<0.50" 4 "<0.75" 5 "<1.00" 6 "<1.15" 7 "<1.25" 8 "<1.50" 9 "<2.00" 10 "<4.00" 11 "4.00+" ) pos(6) row(1) size(small) stack)

graph save "$figures/import_sector_heatmap.gph", replace 


**# Panel B: Initial tariffs

**## Setup

* Import shares	
use "analysis/SIC87_regression_Dataset.dta", clear
keep if year==1978
keep sic ln_imp
replace ln_imp = exp(ln_imp)
collapse (sum) customs  = ln_imp, by(sic)
egen imports_1978 = total(customs)
gen sic_share = customs/imports_1978
keep sic sic_share  cust
gen sic2 = int(sic/100)
collapse (sum)imports = customs sic_share, by(sic2)
merge m:1 sic2 using "processing/sic2_w_desc.dta", keep(1 3) nogen 
save "processing/sic_shares.dta", replace

* Tariffs
use "analysis/SIC87_regression_Dataset.dta", clear
bys sic (year) : gen ord = _n
xtset sic ord
foreach v in year  ave_iv_swiss ave_swiss  ave_ols {
	xtset sic ord
	bys sic (year) : egen `v'_t0 = first(`v')
}
collapse  imports_c2 *_t0, by (sic)
gen sic2 = int(sic/100)

* Merge tariffs and imports		
merge m:1 sic2 using "processing/sic_shares.dta", keep(1 3)
	
* Graph setup
egen first = tag(sic2)
gen timp = imports
gen labpos = 0.0 if first==1
rename sic_share sh_imports

bys sic2 : replace sh_ = . if _n!=1
gsort - sh_
replace sh = sh*100
replace sh = round(sh, 0.05)
gen GLL = sh
tostring sh_imports, replace force
replace sh_ = regexs(1) if regexm(sh_, "((.*)\.[0-9][0-9])")
replace sh = "0" + sh if regexm(sh_, "^\.")
replace sh = sh + "0" if length(sh)==3
replace sh = sh + "%"
replace sh = "" if length(sh) ==3

gen nlab = desc + " (" + sh_imports + ")"
replace nlab = "" if sh==""
replace timp = . if sh==""
gen ord = _n if timp!=.
bys sic2: egen rank_ord = max(ord)
drop ord
rename rank_ord ord

gen ORD = ord+.15

**## Make graph

format ave_ols_t0 % 03.2fc
tw scatter ave_ols_t0 ord , mlc(white) mlp(dash) mc(blue)  msize(medsmall) ///
|| scatter labpos ORD, mlab(nlab) mlabangle(90) m(i)  mlabs(medsmall) mlabc(gs2%80) ///
legend(off) ///
xtitle("") ///
xlabel(none) ///
ytitle("1978 AVE Tariff",size(medium)) ///
xtitle("") ///
caption("(b) Pre-Tokyo 1978 AVE Tariffs by 2-digit SIC Sector", position(6) justification(center) size(medium)  ) ///
ylabel(0[0.05]0.35, labsize(medsmall))

graph save "$figures/initial_tariffs_sic.gph", replace 


	
**# Make combined graph	
	
graph combine   "$figures/import_sector_heatmap.gph" 	"$figures/initial_tariffs_sic.gph", ///
cols(1) /// stack vertically
imargin(medium) /// small gap between panels (top right bottom left)
graphregion(margin(2 2 2 2)) /// remove outer whitespace
plotregion(margin(0 0 0 0)) ///
xsize(4) ysize(6) 
graph export "$figures/fig_1.png", replace	
}
	
****************************
**# Table A1: Summary Stats
****************************
qui {
******************
**# Panel A: SIC 
******************

**## Setup
use "analysis/SIC87_regression_Dataset.dta", clear
keep if year==1979

egen full = rowmiss($tariff_controls $prod_controls $rti_controls ave_ols_ch9 ave_iv_ch9 ln_imp imp_pen ln_exp exp_ship)
keep if full==0


label var ave_iv_ch9 "\$ \Delta \ln (1+AVE_{i}^{IV}) \$"
label var ave_ols_ch9 "\$ \Delta \ln (1+AVE_{i}) \$"
label var exp_ship_ch_9 "\$  \Delta \frac{Exports_{i,t}}{Shipments_{i,t}}   \$"
label var imp_pen_ch_9 	"\$  \Delta \text{Import Penetration}_{i}   \$"
label var ln_exp_ch_9 "\$  \Delta\ln \left(Exports_{i}  \right) \$"
label var ln_imp_ch_9 "\$  \Delta\ln \left(Imports_{i}  \right) \$"
label var ln_skill_pay_diff_ch_9 		"\$ \Delta \ln(\frac{Pay_i^{Non-Prod}}{Pay_i^{Prod}})_{t-1} \$"
label var ln_skill_emp_diff_ch_9 	 	"\$ \Delta \ln(\frac{Emp_i^{Non-Prod}}{Emp_i^{Prod}})_{t-1} \$"
label var ln_wsp_ch_9 			"\$ \Delta \ln(\frac{Wage_i^{Non-Prod}}{Wage_i^{Prod}})_{t-1} \$"

loc vlist  ave_iv_ch9 ave_ols_ch9 exp_ship_ch_9 imp_pen_ch_9 ln_exp_ch_9 ln_imp_ch_9 ///
ln_skill_pay_diff_ch_9 ln_skill_emp_diff_ch_9 ln_wsp_ch_9
 
**## Make table
noi di ""
noi di "**************************************************"
noi di "******************* Table A1 *********************"
noi di "**************************************************"
noi di ""
noi di "***************************"
noi di "******** Panel A **********"
noi di "***************************"
eststo clear
noi tabstat `vlist', s(count mean sd iqr) labelwidth(24) col(stat) 

estpost tabstat  `vlist',  statistics(count mean sd iqr) listwise col(stat) 
esttab using "$results/table_A1_panelA.tex" , cells("count mean(fmt(%9.3f)) sd(fmt(%9.3f)) iqr(fmt(%9.3f))") ///   
collabels("Obs" "Mean" "SD" "IQR") /// 
nostar unstack ///
noobs nonote nonum ///
label replace ///
booktabs ///
substitute(\_ _)


*********************
**# Panel B: Census 
*********************

**## Setup
use "analysis/ipums_industry_data.dta", clear
merge 1:1 ind1990 year using "analysis/ind1990_changes_rout1"
* m=2: inds in ipums emp data that not in industry level data
* 122
* 140
* 301
* 350
* 392

assert _m!=1
keep if _m == 3
drop _merge
drop if inlist(ind1990,210,232,362,390)

loc controls1 
loc controls2 `controls1'	$tariff_controls
di "`controls2'"
loc controls3 `controls2'  	ln_cap_lab_78 skill_emp_share_78 women_ind_share mat_ship_def_78
di "`controls3'"
loc controls4 `controls3'  	$rti_controls
di "`controls4'"
egen full = rowmiss(`controls4')
keep if full==0
drop full

label var ave_iv_ch "\$ \Delta \ln (1+AVE_{i}^{IV}) \$"
label var ave_ols_ch "\$ \Delta \ln (1+AVE_{i}) \$"

label var ln_skill_diff_pay_ch "\$ \Delta \ln(\frac{Pay_i^{Non-Prod,A}}{Pay_i^{Prod,A}})_{t-1} \$"
label var m_ln_skill_diff_pay_ch "\$ \Delta \ln(\frac{Pay_i^{Non-Prod,M}}{Pay_i^{Prod,M}})_{t-1} \$"
label var w_ln_skill_diff_pay_ch "\$ \Delta \ln(\frac{Pay_i^{Non-Prod,W}}{Pay_i^{Prod,W}})_{t-1} \$"

loc vlist  ave_iv_ch ave_ols_ch ln_skill_diff_pay_ch m_ln_skill_diff_pay_ch w_ln_skill_diff_pay_ch

**## Make table
noi di ""
noi di "***************************"
noi di "******** Panel B **********"
noi di "***************************"

eststo clear
noi tabstat `vlist', s(count mean sd iqr) labelwidth(24) col(stat) 

estpost tabstat  `vlist',  statistics(count mean sd iqr) listwise col(stat) 
esttab using "$results/table_A1_panelB.tex" , cells("count mean(fmt(%9.3f)) sd(fmt(%9.3f)) iqr(fmt(%9.3f))") ///   
collabels("Obs" "Mean" "SD" "IQR") ///
nostar unstack ///
noobs nonote nonum ///
label replace ///
booktabs ///
substitute(\_ _)
}
	
******************************
**# Figs A1, 2, C1 & Table C1
******************************
qui{

*************************************
**# Fig A1. Pay inequality densities
************************************* 
		
**## Setup
use "analysis/SIC87_regression_Dataset.dta" if year==1979 | year == 1988, clear
bys sic (year) : gen ord = _n
xtset sic ord

foreach v in  year ln_imp ave_ols ave_iv_swiss ave_swiss  ln_skill_emp_diff ln_skill_pay_diff ln_wsp {  
gen d_`v' = f.`v'-`v'
gen `v'_t0 = `v' if year==1979
}
keep if year==1979  

glob tariff_controls   other_tariff_change  dln_pstar_vw_7279_alt  mfa rho
global prod_controls  ln_cap_lab_78 skill_emp_share_78 women_share 
global rti_controls  ln_invest_def_78 lag_invest_def_ch rti   automation78   mat_ship_def_78
egen full = rowmiss($tariff_controls $prod_controls $rti_controls ave_ols_ch9 ave_iv_ch9  )
keep if full==0

kdensity ln_skill_pay_diff_ch_9, gen(x_pay y_pay)  nodraw
kdensity  ln_skill_emp_diff_ch_9, gen(x_emp y_emp)  nodraw
kdensity ln_wsp_ch_9, gen(x_wsp y_wsp)  nodraw
gen zero = 0 

gen AVE_78 = ave_swiss if year==1979
bys sic : egen ave_78 = max(AVE_78)
drop AVE_78
 
**## Make graph	
format x_* % 03.2fc
format y_* % 01.0fc

tw rarea y_emp zero x_emp , fc(blue%25)  lc(blue) lp(dash)  lw(thin) ///
|| rarea y_wsp zero x_wsp , fc(red%25) lc(red) lp(dash_dot) lw(medthick) ///
||	rarea y_pay zero x_pay , fc(0%20) lc(gs2) lp(solid) lw(thick) ///
legend(order (3 "{&Delta}ln(Pay{sup:Non-Prod} / Pay{sup:Prod})" ///
1 "{&Delta}ln(Emp{sup:Non-Prod} / Emp{sup:Prod})"  ///
2 "{&Delta}ln(Wage{sup:Non-Prod} / Wage{sup:Prod})" ///
) ///
pos(6) ///
row(1) ///
size(medsmall)) ///
xlabel(-1(0.25)1) ///
ytitle("Density")

graph export "$figures/fig_A1.png", width(2000) height(1200) as(png) replace
	
	
 
*****************************************
**# Fig 2. Over vs under liberalization 
*****************************************

**## Setup				
sort ave_swiss_t0
gen yalt = ave_swiss_t0  if mod(_n,5) ==0
gen deviation_col1 = d_ave_swiss - d_ave_ols
gen protected_col1 = deviation_col1<0

gen y1 =  -.1 if _n==1
gen y2 =  -.05 if _n==1
gen x1 =  -.05 if _n==1
gen x2 =  -.1  if _n==1
gen ylab1 = "Over-Liberalized"
gen ylab2 = "Under-Liberalized"
count
loc nobs = r(N)
	
**## Make graph
format d_*  % 03.2fc
format y1 y2 x1 x2 % 03.2fc

tw scatter d_ave_ols d_ave_swiss  , m(i)  mlc(gs2) ///
|| scatter y1 x1 , mlab(ylab1) mlabangle(35) mlabpos(12) m(arrow) msize(vlarge) mc(red%80) msangle(215) mlabc(red%80) mlabsize(large)  yaxis(2) ///
|| scatter y2 x2 , mlab(ylab2) mlabangle(35) mlabpos(6) m(arrow)mc(blue%80) msize(vlarge) msangle(35) mlabc(blue%80) mlabsize(large) yaxis(2) ///
|| scatter d_ave_ols d_ave_swiss if protected_col1==1, mc(blue%20)  mlc(gs2) yaxis(2) ///
|| scatter d_ave_ols d_ave_swiss if protected_col1==0, mc(red%20)   yaxis(2) mlc(gs2)  ///
	, yline(0, lp(dot) lc(red%50)) ///
|| lfitci  d_ave_swiss  d_ave_swiss if d_ave_swiss>-.15, lc(gs2%40) lw(medthick) ///
xtitle("{&Delta}ln(1+AVE{sup:Swiss})") ///
ytitle("{&Delta}ln(1+AVE)", axis(2)) ///
legend(off) 	///
title("")  ///
ylabel(, axis(2)) ///
xlabel(-.2(.05)0, )  ///
ylabel(, labsize(tiny) labcolor(gs2%0) tlc(gs2%0)) ///
yscale( lc(gs2%0) axis(1))

graph export "$figures/fig_2.png", replace width(2000) height(1200)


*****************************************************
**# Fig C1. Densities: over vs under liberalization
*****************************************************

**## Make graph		 
format skill_emp_share % 03.2fc
gen domain = _n/_N
kdensity skill_emp_share if protected_col1 == 1,  gen(x_prot y_prot) nodraw
kdensity skill_emp_share if protected_col1 == 0,  gen(x_nprot y_nprot) nodraw
format y_prot y_nprot x_prot x_nprot % 03.2fc
format y_* % 1.0fc

tw rarea y_prot zero x_prot, m(i) c(l) fcolor(blue%20) lp(dash) lc(blue) ///
|| rarea y_nprot zero x_nprot, m(i) c(l) fcolor(red%20) lc(red) /// 
legend(order(1 "Under-Liberalized" 2 "Over-Liberalized") row(1) pos(6) size(large)) ///
xtitle("Non-Production Share of Workers", size(large)) ///
ytitle("Density", size(large)) ///
xlabel(, labsize(large)) ///
ylabel(, labsize(large))

graph export "$figures/fig_C1.png", replace as(png) width(2000) height(1200)
	

	
******************************************
**# Table C1: Determinants of protection
******************************************

**## Setup	
gen ln_ave_78 = ln(1+ave_78)

label var ln_ave_78 "\$\ln(1+AVE_i)\$"
label var  skill_emp_share_78 "\$ \frac{Emp^{Non-Prod}_i}{Emp_i}\$"
glob tariff_controls other_tariff_change  dln_pstar_vw_7279_alt  mfa rho
global prod_controls ln_cap_lab_78 skill_emp_share_78 women_share 
global rti_controls  ln_invest_def_78 lag_invest_def_ch rti automation78 mat_ship_def_78

eststo clear

eststo:reg 		protected_col1 ln_ave_78, robust
estadd loc ife "N"
estadd loc nobs = e(N)
loc R2  = e(r2_a)
loc R2 :  di %04.3f  `R2'
estadd scalar Rsq = `R2'

eststo:reg 		protected_col1 ln_ave_78	$tariff_controls, robust
estadd loc ife "N"
estadd loc nobs = e(N)
loc R2  = e(r2_a)
loc R2 :  di %04.3f  `R2'
estadd scalar Rsq = `R2'

eststo:reg 		protected_col1 ln_ave_78	$tariff_controls $prod_controls, robust
estadd loc ife "N"
estadd loc nobs = e(N)
loc R2  = e(r2_a)
loc R2 :  di %04.3f  `R2'
estadd scalar Rsq = `R2'

eststo:reg 		protected_col1 ln_ave_78	$tariff_controls $prod_controls $rti_controls, robust
estadd loc ife "N"
estadd loc nobs = e(N)
loc R2  = e(r2_a)
loc R2 :  di %04.3f  `R2'
estadd scalar Rsq = `R2'

**## Make table
noi di ""
noi di "**************************************************"
noi di "****************** Table C1 **********************"
noi di "**************************************************"
noi esttab, se(%9.3f) b(%9.3f) star(* .1  ** .05 *** .01) stats(nobs Rsq)

esttab using "$results/table_C1.tex", ///
noobs ///
nomtitles /// 
label ///
se ///
ar2 ///
b(%04.3f) ///
se(%04.3f) ///
substitute(\_ _) ///
star(* .1 ** .05 *** .01)  ///
nonotes ///
stats( nobs Rsq ,  ///
label("Obs." "Adj. \$R^2\$")) ///
replace		


}


************************************
**# Fig 3: 1st stage & 2nd stage IV 
************************************
qui {
**# Panel A

**## Setup
use "analysis/SIC87_regression_Dataset.dta", clear
keep if year==1979

egen full = rowmiss($tariff_controls $prod_controls $rti_controls ave_ols_ch9 ave_iv_ch9 ln_imp_ch_9  )	 

count
loc nobs = r(N)
format ave* % 03.2fc	
	
**## Make graph
tw scatter ave_ols_ch8 ave_iv_ch9,  mc(blue%80) mlc(white)  ///
|| lfitci ave_ols_ch9 ave_iv_ch9, fc(gs2%40) ,  ///
xtitle("{&Delta}ln(1+AVE{sup:IV})", size(medlarge) margin(b=2)) ///
ytitle("{&Delta}ln(1+AVE)", size(medlarge)) ///
legend(off) ///
caption("(a) First Stage", position(6) justification(center)  size(medlarge)) ///
xlabel(, labsize(medlarge)) ///
ylabel(, labsize(medlarge)) ///

graph save "$figures/ave_vs_ave_iv_c2.gph", replace

	
**# Panel B	

**## Setup
reg ave_ols_ch9 ave_iv_ch9 
predict ann_dave_hat
format ann* % 03.2fc
count
loc nobs = r(N)
format ann* % 03.2fc
	
**## Make graph
tw scatter ln_imp_ch_9 ann_dave_hat  if !inlist(sic, 2024, 2075)  , mc(blue%80) mlc(white) /// 
|| lfitci ln_imp_ch_9  ann_dave_hat  ,  fc(gs2%40) clcolor(green%60)   ///
xtitle("Predicted {&Delta}ln(1+AVE)",  size(medlarge) margin(b=2)) ///
ytitle("{&Delta}ln(Imports)", size(medlarge)) ///
caption("(b) Second Stage", position(6) justification(center) size(medlarge)  ) ///
legend(off) ///
xlabel(, labsize(medlarge)) ///
ylabel(, labsize(medlarge)) ///

graph save "$figures/dln_imports_vs_ave_hat.gph", replace
	

**# Make combined graph	
graph combine   "$figures/ave_vs_ave_iv_c2.gph" 	 "$figures/dln_imports_vs_ave_hat.gph" 
graph export "$figures/fig_3.png", replace width(2000) height(1200)
}		
		

*******************************************
**# Fig A5: SIC v IPUMS dataset comparison  
*******************************************
qui {
**# Setup
use "analysis/SIC87_regression_Dataset.dta", clear

keep if year==1979

egen full = rowmiss($tariff_controls $prod_controls $rti_controls ave_ols_ch9 ave_iv_ch9 ln_imp    )
keep if full==0

sum ave_ols_ch9,d
glob ave_s  `=`r(p75)'-`r(p25)'' 
sum ln_skill_pay_diff_ch_9, d
glob skill_pay_s  `=`r(p75)'-`r(p25)'' 
gen obs = _n

keep obs ave_ols_ch9 ln_skill_pay_diff_ch_9 
tempfile somedata
save `somedata', replace

use "analysis/ipums_industry_data.dta", clear
merge 1:1 ind1990 year using "analysis/ind1990_changes_rout1"
* m=2: inds in ipums emp data that not in industry level data
* 122
* 140
* 301
* 350
* 392
assert _m!=1
keep if _m == 3
drop _merge
drop if inlist(ind1990,210,232,362,390)

sum ave_ols_ch,d
glob ave_i  `=`r(p75)'-`r(p25)'' 
sum ln_skill_diff_pay_ch,d
glob skill_pay_i  `=`r(p75)'-`r(p25)'' 

gen obs = _n
merge 1:1 obs using `somedata', nogen

sum ln_skill_pay_diff_ch, d
sum ln_skill_diff_pay_ch, d

gen xvals_skill =-.82 +  _n*(1/2)/30  
replace xvals_skill = . if !inrange(xvals_skill, -.82,.9)

kdensity ln_skill_diff_pay_ch, nograph at(xvals_skill) gen(pay_x_i pay_y_i)
kdensity ln_skill_pay_diff_ch, nograph at(xvals_skill) gen(pay_x_n pay_y_n)

gen xvals_ave = -0.145 + _n*(1/2)/200
replace xvals_ave = . if !inrange(xvals_ave, -0.10, 0.08)

kdensity ave_ols_ch, nograph at(xvals_ave) gen(ave_x_i ave_y_i)
kdensity ave_ols_ch9 , nograph at(xvals_ave) gen(ave_x_n ave_y_n)

format xvals_ave xvals_skill %03.2fc

gen zero = 0

 
**# Make graph 
* Skill pay
tw rarea pay_y_n zero xvals_skill 	, fc(gs2%20) lc(gs2%40) lp(dash) ///
|| rarea pay_y_i zero xvals_skill 	, fc(purple%20) lc(purple%40)   ///
xtitle("{&Delta}ln(Pay{sup:Non-Prod}/Pay{sup:Prod})", size(large)) ///
xlabel(, labsize(large)) ylabel(, labsize(large)) ///
legend(order( 1 "SIC"  2 "IPUMS") row(1) pos(6))  ///
ytitle(Density, size(large))
graph save "$figures/Dataset_comparison_skillpay.gph", replace

* AVE
tw rarea ave_y_n zero xvals_ave 	, fc(purple%20) lc(purple%40)   ///
|| rarea ave_y_i zero xvals_ave 	, fc(gs2%20) lc(gs2%40) lp(dash) ///
xtitle("{&Delta}ln(1+AVE)", size(large)) ///
xlabel(, labsize(large)) ylabel(, labsize(large)) ///
legend(order( 1 "SIC"  2 "IPUMS") row(1) pos(6) size(large)) ///
ytitle(Density, size(large))
graph save "$figures/Dataset_comparison_ave.gph", replace

* Make graph
grc1leg "$figures/Dataset_comparison_skillpay.gph" ///
"$figures/Dataset_comparison_ave.gph", ///
row(2)
graph export "$figures/fig_A5.png", replace 
}	

log close
		